Insights into particle dispersion and damage mechanisms in functionally graded metal matrix composites with random microstructure-based finite element model

This study investigates the impact of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathrm {Al_2O_3}$$\end{document}Al2O3 particle volume fraction and distribution on the deformation and damage of particle-reinforced metal matrix composites, particularly in the context of functionally graded metal matrix composites. In this study, a two-dimensional nonlinear random microstructure-based finite element modeling approach implemented in ABAQUS/Explicit with a Python-generated script to analyze the deformation and damage mechanisms in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathrm{AA6061\mbox{-}T6/Al_2O_{3}}$$\end{document}AA6061-T6/Al2O3 composites. The plastic deformation and ductile cracking of the matrix are captured using the Gurson–Tvergaard–Needleman model, whereas particle fracture is modelled using the Johnson–Holmquist II model. Matrix-particle interface decohesion is simulated using the surface-based cohesive zone method. The findings reveal that functionally graded metal matrix composites exhibit higher hardness values (\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textrm{HRB}$$\end{document}HRB) than traditional metal matrix composites. The results highlight the importance of functionally graded metal matrix composites. Functionally graded metal matrix composites with a Gaussian distribution and a particle volume fraction of 10% achieve \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textrm{HRB}$$\end{document}HRB values comparable to particle-reinforced metal matrix composites with a particle volume fraction of 20%, with only a 2% difference in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textrm{HRB}$$\end{document}HRB. Thus, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textrm{HRB}$$\end{document}HRB can be improved significantly by employing a low particle volume fraction and incorporating a Gaussian distribution across the material thickness. Furthermore, functionally graded metal matrix composites with a Gaussian distribution exhibit higher \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textrm{HRB}$$\end{document}HRB values and better agreement with experimental distribution functions when compared to those with a power-law distribution.

Functionally graded metal matrix composites ( FGMMCs ) belong to the category of metal matrix composites ( MMCs ) and are characterized by a continuous variation in the volume fraction ( VF r ) of the reinforcement along a specific direction of the matrix alloy 6 .FGMMCs offer a gradual or continuous transition in engineering properties at the macroscopic scale, allowing for the combination of desirable properties without the presence of mechanically weak interfaces, which is a limitation often encountered in surface coating techniques 7 .
A literature review revealed the existence of multiple particle distribution curves within FGMMCs ; examples are: [8][9][10][11][12][13][14][15][16] , as depicted in Fig. 1.These empirical curves serve as a valuable resource for theoretical investigations, enabling the development of a controllable function that accurately represents particle distribution within the composite material.Utilizing such a controllable function facilitates a systematic approach for manipulating and analyzing the particle distribution.By fitting the experimental distribution curve to an approximated function, it is possible to obtain a mathematical model that captures the essential characteristics of particle distribution in FGMMCs .This model can be integrated into theoretical studies to enable a deeper understanding of the behavior of FGMMCs , which can serve as a powerful tool for predicting and optimizing the performance of FGMMCs in practical applications.
To assess the mechanical properties of FGMs , the indentation method is often favored owing to its simplic- ity and minimal specimen preparation requirements.This technique can be easily implemented multiple times on both both small-scale materials and miniaturized structures using a single specimen and a suitable load and indenter tip geometry.The use of numerical modeling to simulate the indentation process of FGMMCs can provide insight into the deformation mechanisms, stress, and strain fields within the metal matrix surrounding the reinforcement particles, as well as the interaction between the particles and the matrix 17 .
Reference 18 conducted a study wherein graded Al7075/SiC composites were fabricated using a centrifugal casting technique with the aim of enhancing the mechanical properties and wear resistance for automotive purposes.The fabrication process involved utilizing two distinct weight fractions of particles (6.5% and 9.5%), while maintaining a constant mold rotational speed of 1300 rpm.The findings of the study indicated an improvement in the mechanical properties and wear resistance of the composites, especially in the outer region, with an increase in the weight fraction of particles owing to the influence of centrifugal forces.
In the study conducted by 19 , the distribution of SiC particles in aluminum-based FGMMCs used in brake rotor discs was investigated.The average particle size used in this study was 23 µm .The results demonstrated that the maximum hardness achieved after heat treatment of the Al(356)-SiC and Al(2124)-SiC FGMMCs at the outer periphery was 155 BHN and 145 BHN, respectively.Reference 20 conducted a study on Al-SiC FGMMCs , where they varied the volume percentage of reinforcement.They found that, as the percentage of reinforcement increased, the density decreased.However, they observed that the hardness and wear resistance increased from the core to the cast surface.At low sliding speeds, they noted that microcracking and abrasive wear were dominant factors.
Reference 21 conducted an investigation into the effects of particle segregation ratio and particle distribution on the VF r and particle size ( PS ) of reinforcement in aluminum based FGMMCs with in-situ primary Si/Mg 2 Si particles.The study involved varying parameters, such as the mold temperature and high pouring temperature.The fabricated FGMMCs showed a graded dispersion of primary Si and Mg 2 Si particles in the inner region, with diameters ranging from 70 to 30 µm and 30 to 18 µm , respectively.Additionally, the volume fractions of Mg 2 Si and primary Si exhibited a gradient distribution ranging from VF r of 14.8 to 27.7%.The highest hardness values were observed in the middle portion of the inner layer of the fabricated FGMMCs , ranging from Rockwell hardness ( HRB ) 72.0 to 75.0.Conversely, the outer layer exhibited the lowest hardness value of approximately HRB 56.0 63.0.
Reference 22 enhanced the resistance to contact damage in a graded material system comprising silicon nitride ( Si 3 N 4 ) ceramic (higher modulus) and oxynitride glass (lower modulus) by employing a unidirectional gradient of the elastic modulus from the contact surface to the interior.They achieved a 30% reduction in the maximum tensile stresses outside the Hertzian contact circle compared with the monolithic material ( Si 3 N 4 ) of the graded system.
Reference 23 investigated the effectiveness of two-scale modeling for analyzing the mechanical properties of FGMs and evaluated aluminum as the supporting matrix and silicon carbide as particle inclusions.The mechani- cal properties of the material were derived at the macroscale level using a full model, and at the microscale level using a representative volume element in two-scale modeling.The results demonstrated the usefulness and reliability of the two-scale model in reducing the numerical computational time.As the VF r of the inclusion increases, the deviation in the stress values decreases, providing evidence supporting the theory that FGMs offer the advantage of smoother stress distributions.
Reference 17 conducted a study on the indentation behavior of FGMMCs , revealing that the number of layers, compositional gradient exponent, and random particle dispersion have a significant impact on the properties of material.Increasing the number of layers resulted in noticeable increases in indentation depths, whereas increasing the compositional gradient exponent led to higher mean residual stresses.Conversely, at a specific layer number, an increase in the compositional gradient exponent decreased the mean residual stresses and strains owing to an increase in the ceramic VF r .Random particle dispersion influences the central indentation depth and deformed surface profiles, resulting in non-uniform levels and distributions of residual stress and strain.
To the best of our knowledge, no computational FE model has been developed thus far to quantitatively analyze the evolution of multiple damage mechanisms during the indentation of FGMMCs .The development of such a model is crucial for predicting material behavior under various conditions by systematically analyzing failure mechanisms and simulating responses to stress and strain over time.These models are essential for optimizing the material performance because they enable the tailoring of microstructures to enhance durability and resilience, thus reducing the need for extensive physical testing.They manage complex geometries and loading conditions, providing a virtual testing environment that identifies potential failure points and predicts the realworld performance.Integrating finite element damage models into research and development not only improves cost efficiency and innovation but also deepens material understanding and facilitates advancements in material engineering across various industries.The objective of this study is to address the gap in stress-state-dependent quantitative damage growth analysis in FGMMCs by investigating the indentation of AA6061-T6/Al 2 O 3 FGMMCs .To achieve this, a 2D microstructure-based FE model is employed to closely approximate the actual microstructure and consider all three potential failure modes.The hardness assessment procedure outlined in the ASTM E18-22 standard 24 is incorporated into finite element analysis ( FEA ).To capture the damage mechanisms in the AA6061-T6 matrix, the Gurson-Tvergaard-Needleman ( GTN ) model is employed.This model represents the behavior of the matrix considering void growth and coalescence.Simultaneously, the Johnson-Holmquist II ( JH2 ) model is integrated into the FE model to represent Al 2 O 3 particle cracking.The matrix-particle interfacial behavior is further modeled using the cohesive zone method ( CZM ).The analysis is conducted using ABAQUS/ Explicit FE software, complemented by Python for scripting and preprocessing 25,26 .

Problem formulation
Figure 2 depicts the configuration of an AA6061-T6/Al 2 O 3 FGMMCs with a spherical indenter.The composite structure is simulated as a rectangle with a width ( W ) of 5 mm and height ( H ) of 2.5 mm .In accordance with the ASTM E18-22 standard 24 , the diameter of indenter ( D ind ) and applied force ( F ) were determined as 1.558 mm and 980 N, respectively.The HRB was computed using the following equation where h represents the permanent penetration in millimeters (mm).

Material constitutive models
This study focuses on investigating a composite material consisting of AA6061-T6 alloy, which is chosen as the base matrix because of its extensive range of engineering and structural applications 7 .The composite is reinforced with Al 2 O 3 particulates, which are the most commonly used reinforcement for PRMMCs 27 with volume fractions of 10% and 20% and an average size of 100 µm .In general, composite materials can fail because of one of the three major failure modes.These modes include plastic deformation and ductile cracking of the matrix, fracture of the reinforcement particles, and decohesion at the matrix-particle interface.The following subsections provide a detailed description of each failure mode.

Plastic deformation and ductile cracking of the matrix
A three-stage mechanism involving void nucleation, growth, and coalescence is responsible for the ductile damage observed in the metals.These voids can be initiated by inclusions or microcracks.Subsequently, these voids undergo enlargement owing to the accumulation of plastic strain.The GTN model 28 is commonly employed for the theoretical analysis of this phenomenon in porous metals.The GTN model is a modified version of the Gurson model 29 and its yield function is mathematically expressed as where the non-dilatational strain energy is denoted by ϕ , and constants q 1 , q 2 , and q 3 are introduced by Tvergaard 30 to account for void interactions.Von Mises and flow stresses of the intact material are represented by σ q and σ y , respectively.
To accurately capture the rapid reduction in stress-carrying capacity resulting from void coalescence, Tvergaard-Needleman 28 introduced the parameter f * , known as the effective porosity.This parameter acts as a modeling tool and is mathematically defined as follows: The GTN model incorporates several parameters to describe void behavior.The critical void volume fraction ( VVF ) at the beginning of coalescence is denoted by f c and f * u = 1 q 1 represents the VVF at which the material loses its stress-bearing capacity.Meanwhile, f F represents the VVF at which the material experiences complete failure and controls the element deletion process.The increase in the VVF owing to void nucleation and growth is considered.The function describing the effective porosity can be expressed as follows: The rate of void growth ( df g ) can be mathematically represented as a function of plastic volume change when the material is considered plastically incompressible.
where the trace of the plastic strain-rate tensor is expressed by dε P ii .The nucleation of voids is strongly influenced by plastic strain, particularly under hydrostatic tension 31 .This phenomenon can be mathematically represented as follows Variable p represents the hydrostatic stress, and f N represents the volume fraction of nucleated voids.Addition- ally, ε N signifies the mean equivalent plastic strain for void nucleation and s N represents the standard deviation of the void distribution.The macroscopic plastic work rate can be equivalent to the rate of matrix plastic dissipation, which in turn determines the rate of equivalent plastic strain dε p according to the following equation Table 1 presents the elastoplastic and GTN model constants employed for the AA6061-T6 matrix derived from previous studies.

Fracture of reinforcement particles
The JH2 model constitutive model was initially employed to simulate the mechanical response of materials exhibiting brittle fractures, particularly ceramic materials.Based on the foundational principles outlined in 32 , the JH2 model incorporates the mechanisms of softening and pressure-dependent strength, material damage, and fracture; it also accounts for significant residual strength after fracture, bulking, and sensitivity to the loading rate.The JH2 model characterizes the behavior of materials in their damaged state by considering three distinct states: intact, damaged, and fractured.In this context, the model employs an analytical function to express the normalized equivalent stress in the damaged state.The generic form of this function is as follows , and σ HEL are defined as follows: σ * i denotes the normalized intact equivalent stress; σ * f represents the normalized fracture stress; D is the damage factor, which varies between 0 and 1.0; and σ HEL denotes the equivalent stress at the Hugoniot Elastic Limit ( HEL ).The critical point signifies the net compressive stress, which considers both hydrostatic pressure and deviatoric stress components.This is the point at which a onedimensional shock wave with uniaxial strain exceeds the elastic limit of the material 33 .Lastly, σ denotes the actual equivalent stress.
The normalized intact strength ( σ * i ) can be determined using the following equation: where A, C, and N are the material constants.P * represents the normalized pressure defined as the ratio of the actual hydrostatic pressure (P) to the hydrostatic pressure at the HEL ( P HEL ).T * represents the normalized maximum tensile hydrostatic pressure, which is defined as the ratio of the maximum tensile hydrostatic pressure that the material can withstand (T), to P HEL .ε * represents the dimensionless strain rate, defined as the ratio of the actual equivalent strain rate ( ε ) to the reference strain rate ( ε0 = 1.0 s −1 ).Similarly, the equation for the normalized fracture strength ( σ * f ) can be expressed as where B, C, and M are the material constants.SFMAX represents the ultimate value of the normalized fracture strength ( σ * f ), providing additional flexibility in the definition of the fracture strength.According to the JH2 model, the softening process in brittle materials can be described by Eq. ( 9), which allows for the gradual softening of the material as the plastic strain increases.The softening process continues until the material is fully damaged ( D p = 1 ).The expression describing the cumulative damage resulting from fracture is expressed as follows ( 9) where �ε P represents the accumulated plastic strain during the integration cycle.Function ε P f = f (P) represents the plastic strain necessary for fracture under a constant pressure P. The parameters d 1 and d 2 correspond to the damage factors associated with ε P f .When a material reaches a certain threshold of plastic deformation or damage, it enters a failure state described by fluid-like behavior 33 .In this state, the material loses its strength and cannot withstand the stress.Both the hydrostatic pressure and deviatoric stress become zero.The relationship between hydrostatic pressure P and volumetric strain µ is described by a polynomial equation of state ( EOS ), which consists of two distinct stages: an elastic stage and a plastic damage stage.The mathematical expressions for these stages are as follows Here, K 1 , K 2 , and K 3 are constants and µ is defined as the ratio of the current density ( ρ ) to the initial density ( ρ 0 ).For tensile pressure ( µ < 0 ), Eq. ( 13) is replaced by P = K 1 • µ .The incremental pressure ( ∆P ) is added when the material fractures owing to bulking energy.
The reduction in incremental internal elastic energy is transformed into potential internal energy via an incremental increase in hydrostatic pressure ∆P .As the fracture progresses, the shear and deviatoric stresses diminish owing to the decrease in the equivalent plastic flow stress σ .The elastic internal energy related to these shear and deviatoric stresses can be expressed mathematically as where G denotes the rigidity modulus.The incremental energy loss is given by The quantities U D(t) and U D(t+�t) are determined using Eq. ( 14).The change in energy, ∆U , is primarily con- verted into incremental fracturing energy, ∆F .An approximate expression for this energy conversion is given by where β is a fraction satisfying (0 ≤ β ≤ 1) , which represents the extent of energy transformation.The JH2 model constants used for the Al 2 O 3 particles are summarized in Table 1.

Cohesive zone modeling (CZM) for matrix and particle interfaces
The failure process of an interface surface involves three fundamental components: initiation of damage, progression of damage, and overall debonding resulting from substantial damage.The cohesive zone method ( CZM ) is utilized to describe the potential debonding occurring at the interface between the particles and matrix.This study characterizes the CZM through a bilinear relationship between traction ( T ) and separation ( ).The initiation of damage is governed by the quadratic stress failure criterion.This criterion postulates that damage is initiated when a quadratic interaction function reaches a unit value 38 .Mathematically, this criterion can be expressed as where T 0 n , T 0 s , and T 0 t represent the peak values of the nominal stresses during deformation along specific directions: normal to the interface, the first shear direction, and the second shear direction, respectively.
In the traction-separation model, it is assumed that the scalar damage variable D changes from 0 to 1 as the material undergoes further loading after damage initiation.The damage variable D exerts the following influence on the stress components: The stress components σ * n , σ * s , and σ * t are obtained by utilizing the elastic traction-separation behavior for the current strain before damage initiation.To characterize the progression of damage within the interface under combined normal and shear deformations, it is helpful to introduce an effective displacement m .The effective displacement is defined as follows ( 12) where 0 m and f m characterize the effective separations at the initiation of damage and complete failure, respectively, and max m represents the peak value of the effective displacement attained during the loading history.The fracture energy, G c , which is a measure of the area enclosed by the traction-separation displacement curve, can be obtained as follows The CZM constants used for the AA6061-T6 and Al 2 O 3 interfaces in this study are listed in Table 2.

Computational model
The FE model and its validation are described in the following subsection.To investigate the damage mechanisms related to the failure modes of AA6061-T6/Al 2 O 3 composites, a Python-generated script is used.This script enables the implementation of indentation-hardness test specifications to assess the resulting microstructure.Using this script, it is possible to generate particles with different distributions and volume fractions throughout the FGMMCs structure.
Mathematical equations representing the different failure modes outlined in this study have been compiled and integrated into the ABAQUS FE software.The analysis process is summarized in the flowchart presented in Fig. 3, providing an overview of the steps involved in the current study.

Finite element model
The analysis considered various factors such as geometric considerations, material nonlinearities, and nonlinear deformations caused by the contact between a rigid indenter and a metal matrix surface.The analysis employed a 2D FE model with 4-node bilinear plane strain quadrilateral elements ( CPE4R ).To capture the damage behav- ior and failure modes precisely, a refined mesh with an element size of 0.01 mm is employed throughout the microstructure.This resulted in an average of 146,700 elements and 156,000 nodes per microstructure (Fig. 4).In accordance with the findings of a previous investigation 41 , the optimal particle shape utilized in this study is circular with a diameter (d) of 100 µm To optimize computational efficiency, the FE model employs a perfectly rigid indenter, and a fixed lower boundary is chosen for the structure to maintain stability.A Coulomb friction coefficient of 0.1 42,43 is applied at both the indenter-matrix and particle-matrix interfaces to model frictional interactions.For quasi-static analysis in the ABAQUS software, the force was gradually applied in a smooth step 44 , as depicted in Fig. 5.

Validation of the model
The accuracy of the FE model is validated by comparing its predicted HRB values with those obtained experi- mentally for the corresponding PRMMCs reported by 45 .The measurement procedure defined in ASTM E18-22 24 is employed in the experimental study.This method initially involves applying a minor load to establish a reference position, followed by the application of a major load for a specified time interval.Subsequently, the major load is removed and the difference in the penetration depth of the indenter is measured and used to calculate the HRB value.
For both volume fractions of 10% and 20%, at least ten random base models were considered in the experimental study, as the Al 2 O 3 particles were randomly dispersed within the structure (refer to Fig. 6).The average HRB value was used for the comparison.
Table 3 presents a comparison between the random distribution of the proposed model and experimental data.The results indicate that the error ranges from 0.635 to 1.811% when the present random distribution is compared with the experimental data.This low percentage error indicates that the HRB values obtained from the FE model for volume fractions of 0%, 10%, and 20% are very close to the corresponding experimental values.The comparison demonstrates good agreement between the predicted and experimental results.This minimal error percentage error suggests that the HRB values derived from the FE model for volume fractions of 0%, 10%, and 20% closely align with the corresponding experimental values.The comparison demonstrates strong agreement www.nature.com/scientificreports/ between the predicted and experimental results.For further details about the model and its validation, please refer to the previous work 41 .

Results and discussions
The main objective of this section is to explore the effect of Al 2 O 3 particle distribution on the mechanical char- acteristics and damage mechanisms of the AA6061-T6/Al 2 O 3 composite.To achieve this objective, this study analyzes the effects of the Al 2 O 3 particle distribution across the FGMMCs structure on the mechanical properties of the composite.An analysis is conducted to investigate the indentation behavior of structures composed of FGMMCs under the influence of a spherical indenter.The structure under examination is modelled to exhibit a variation in the VF r of the Al 2 O 3 particles throughout its thickness, as depicted in Fig. 2.
FGMMCs are modelled using the linear rule of mixtures, which is a commonly employed approach 17 .This modelling technique assumes a linear relationship between the volume fractions of the Al 2 O 3 reinforcement and AA6061-T6 matrix material.The sum of the volume fractions of the reinforcement and matrix material is equal to one and is expressed as  www.nature.com/scientificreports/V rein and V mat represent the volume fractions of Al 2 O 3 and AA6061-T6 matrix materials, respectively.In the investigated FGMMCs structures (Fig. 2), the Al 2 O 3 particles are distributed randomly throughout the matrix.This irregularity in the arrangement of the reinforcement particles reflects the actual conditions observed in the FGMMCs.
The investigation focused on studying the variation in VF r of Al 2 O 3 particles across the thickness of FGMMCs using two approaches.The first approach is the power law approach 46 , whereas the second approach utilized a Gaussian distribution 47 .
The first approach (power law) is commonly employed in the theoretical modeling of FGMMCs structures 17,48,49 .In this approach, the VF r of Al 2 O 3 particles at any position y across the thickness is assumed to be determined by the equation Here, n represents the non-negative compositional gradient exponent, which determines the variation profile of the Al 2 O 3 particles.Different values of n result in distinct particle distribution profiles and overall particle volume fractions.Position y indicates the location within the thickness of the FGMMCs structure measured from its bottom surface.H represents the total thickness of the FGMMCs structure.
To enable meaningful comparisons between PRMMCs and FGMMCs , two total particle volume fractions are considered, 10% and 20%.To determine the volume fractions, the compositional gradient exponent n is   www.nature.com/scientificreports/estimated to be 2.65 for a total VF r of 10%, whereas a value of 0.65 is estimated for n when the total VF r is 20%.
Figure 7 shows the distribution curves obtained using the estimated values.Figure 8 shows the samples representing FGMMCs resulting from the application of the power law approach.A sample with a total VF r of 10% ( n = 2.65 ) is depicted in Fig. 8a, whereas a sample with a total VF r of 20% ( n = 0.65 ) is illustrated in Fig. 8b.
The second approach (Gaussian distribution) used to vary the Al 2 O 3 particle distribution within the thickness of the FGMMCs structure is based on experimental work, specifically the study conducted by [11][12][13]16 (see Fig. 1). In his approach, the VF r distribution at any position y is determined fitting the experimental distribution curves.The best fit is achieved using a Gaussian distribution function 47 , which is represented by the following formula: The constants A 0 , x 0 , and σ G determine the shape of the Gaussian curve, and consequently, the particle volume fraction.These constants are estimated using trial methods to obtain specific volume fractions, these constants are estimated through trial methods.For a total VF r of 10%, the estimated values of A 0 , x 0 , and σ G are 0.365, 2.215, and 0.295, respectively.Similarly, for a total VF r of 20%, the constants are estimated to be A 0 = 0.365 , x 0 = 1.95 , and σ G = 0.65 .The distribution curves for VF r values of both 10% and 20% using a Gaussian distribu- tion are shown in Fig. 9.
Representative samples of FGMMCs structure with these VF r variations obtained using the Gaussian approach are shown in Fig. 10. Figure 10a depicts a structure with a total VF r of 10%, whereas Fig. 10b shows a structure with a total VF r of 20%.

HRB results
Nonlinear FE analyses are performed to explore the influence of the power law and Gaussian distributions in FGMMCs on the indentation behavior of AA6061-T6/Al 2 O 3 FGMMCs .Specifically, two total volume fractions of Al 2 O 3 particles are considered: 10% and 20%.This is achieved by varying the compositional gradient expo- nents ( n = 0.65 and 2.65 ) for the power law distribution and by adjusting three parameters ( A 0 , x 0 , and σ G ) for the Gaussian distribution.
Ten samples with random reinforcement distributions are analyzed for each VF r and distribution, and the calculated indentation depths are used to determine the HRB values.Figure 11 presents the HRB results for the  pure AA6061-T6 matrix as a reference value and the average HRB values for the AA6061-T6/Al 2 O 3 composites with volume fractions of 10% and 20% Al 2 O 3 particles.In addition to the HRB of the pure matrix (0% VF r ), the HRB results are shown three times in Fig. 11 for both the 10% and 20% particle VF r .The first set of results corresponds to the random distribution of Al 2 O 3 particles throughout the thickness without any variation in par- ticle concentration, labeled as PRMMCs .The other two sets of HRB results show the effects of particle variation through the thickness using the power law and Gaussian distribution, labeled FGM-1 and FGM-2, respectively.www.nature.com/scientificreports/According to the results presented in Fig. 11, the HRB results exhibited a consistent pattern.Specifically, when examining volume fractions of 10% and 20%, it is apparent that the HRB values of PRMMCs are greater than those of the pure AA6061-T6 matrix (0% VF r ).Additionally, it can be observed that the HRB values of the FGMMCs using the power law distribution (labeled FGM-1) are higher than those achieved with the PRMMCs , and the HRB values of FGMMCs obtained through the Gaussian distribution (labeled FGM-2) are higher than those obtained through the power law distribution.
Table 4 shows a comparison between the HRB results and their enhancement over the pure AA6061-T6 matrix for all Al 2 O 3 particle volume fractions and variations through the structure thickness.It is shown that the HRB is enhanced for a 10% VF r over the pure matrix by 11.7%, 15.6%, and 17.8% for PRMMCs , FGM-1, and FGM- 2, respectively.For 20% VF r , the enhancements are 19.8%, 25.2, and 25.7 for PRMMCs , FGM-1, and FGM-2, respectively.It is noted from the results in Table 4 that the hardness values in the case of the Gaussian distribution are greater than in the case of the power law distribution as a result of the convergence of the particles in the part closest to the surface of the samples, which leads to an increase in the resistance of the samples to penetration due to the increase in the hardest material in the upper part of the samples, this difference in results appears as a result of the difference in the distribution according to the shape of the functions, which is clearly shown in Figs.7 and 9, as well as the distribution in samples representing the distributions in Figs. 8 and 10.
The results highlight the importance of FGMMCs , because the HRB values obtained for FGMMCs with a Gaussian distribution and a particle VF r of 10% are comparable to the HRB values obtained for PRMMCs with a particle VF r of 20%.The difference between the HRB results is approximately 2%.This means that with a low particle VF r and changing the particle variation across the material thickness, that is, by a Gaussian distribution, HRB can be enhanced significantly.
As particle VF r increases, it becomes evident that the HRB values also increase.Notably, the results obtained for both distributions of FGMMCs with a particle VF r of 20% exhibited minimal differences.This observation indicates that it is challenging to achieve a significant difference in the HRB values through particle distribution for FGMMCs with higher particle volume fractions such as 20%.
The load-displacement and surface profiles curves of representative specimens are shown in Fig. 12.The indentation depth decreased significantly as VF r increased.Notably, distinct variations in the results are evident for all representative samples with a VF r of 10%, as illustrated in Fig. 12a,c.Conversely, for a VF r of 20%, the curves for the two representative samples of FGMMCs are nearly identical, as depicted in Fig. 12b,d.

Residual von-Mises stress and plastic strain
Figure 13 illustrates the effect of VF r on the distributions of residual von-Mises stress and effective plastic strain within the FGMMCs structure.As VF r decreases, the levels of residual stress and strain in the representative sam- ples decrease.It is worth noting that the regions of the effective plastic strain are localized beneath the indenter and propagate at an angle of nearly 45 • .
The stress concentration is influenced by the relative distances between the particles and indenter as well as between the particles themselves.Consequently, as VF r increases, FGMMCs becomes more susceptible to stress concentration.This is because the permissible area for particle movement is reduced with higher particle VF r .Moreover, the particles located beneath the indenter experienced the most significant effects, as shown in Fig. 13b,d.
The range of von-Mises stress varies from 1010 MPa to 1565 MPa for a particle VF r of 10% ( n = 2.65 ), and for a particle VF r of 20% ( n = 0.65 ), with a range of 1282 MPa to 3215 MPa , as shown in Table 5.The effective plastic strain ranges from 1.566 to 2.41 for a VF r of 10% ( n = 2.65 ), and 1.507 to 2.524 for VF r of 20% ( n = 0.65 ) (see Table 5).
Figure 14 presents the effects of the Gaussian distribution on the distributions of the residual von-Mises stress and effective plastic strain within the FGMMCs structure for both 10% and 20% volume fractions.Similar to previous observations, localized stress concentration areas beneath the indenter and the propagation angle of the effective plastic strain are evident in these results (see Fig. 14b,d).
By analyzing the distribution of von-Mises stress for the representative samples (Fig. 14a,b), it is observed that the stress levels are higher for 10% VF r .This can be attributed to the chosen Gaussian distribution, which results in a higher concentration of particles in the region directly below the indenter compared with the distribution with a 20% VF r .This can be observed by noting the positions of the peak values in the curves, as depicted in Fig. 9.
The range of von-Mises stress for FGMMCs with a Gaussian distribution is 1401 MPa to 2446 MPa for a particle VF r of 10% and 1196 MPa to 2228 MPa for a particle VF r of 20%, as shown in Table 6.The effective plastic strain ranged from 1.615 to 2.865 for the 10% particle VF r and 1.331 to 2.474 for the 20% particle VF r (see Table 6).These values indicate the extent of plastic deformation experienced by the matrix material in the FGMMCs structures with Gaussian distributions.A comparison of the von-Mises stress results between FGMMCs and PRMMCs showed that there is a signifi- cant increase in the von-Mises stress values in the case of FGMMCs .For instance, when comparing the average von-Mises stress values in both cases, it is found that with 10% VF r , the von-Mises stress values increased by 64% and 117% for FGMMCs with power law and Gaussian distributions, respectively.Similarly, with 20% VF r , the average value of the von-Mises stress increased by 81% and 57% for both the power law and Gaussian distributions, respectively.
Furthermore, an increase in the average values of effective plastic strain is also observed.When considering a 10% VF r of particles, the effective plastic strain increased by 84% and 90% for the power law and Gaussian distri- bution, respectively.Similarly, with 20% VF r , the effective plastic strain for power law and Gaussian distribution increased by 35% and 37%, which indicates the tendency to damage for FGMMCs over PRMMCs.

Damage analysis of FGMMCs
This subsection presents an analysis of the existing and potential damage locations within the AA6061-T6/Al 2 O 3 FGMMCs structure for three potential failure modes: (i) matrix damage, (ii) particle fracture, and (iii) interfacial decohesion.The results of the damage parameters associated with each failure mode are analyzed.
(i) Matrix damage: For VF r of 10%, no damage is observed in any sample for both FGMMCs distribution approaches.Fig. 15a,c shows the representative samples for this VF r , which help predict the failure mechanisms in the AA6061-T6 matrix.As shown in Table 5, the range of the void volume fraction ( VVF ) for a VF r of 10% ( n =2.65) is between 0.00121 and 0.00398, with an average value of 0.00216.Considering the results of VVF for the Gaussian FGMMCs distribution, the range of VVF is between 0.00115 and 0.01254, with an average value of 0.00299 for a VF r of 10% (see Table 6).
Based on the results, it is observed that the matrix damage mode occurred when VF r is 20% for both the power law and Gaussian distributions of the FGMMCs .Out of the ten samples studied for each distribution, only one sample in each case exhibited clear damage.These two samples ( S 2 and S 7 ) are shown in Fig. 15b,d.For a VF r of 20% ( n =0.65), the range of the VVF is between 0.00229 and 0.06, with an average value of 0.012353 (Table 5).These findings indicated that a higher VF r leads to a higher possibility of void nucleation and growth.For a VF r of 20% Gaussian FGMMCs distribution, the range is between 0.00162 and 0.06, with an average value of 0.00913 (see Table 6).The higher possibility of matrix damage in these cases is attributed to the increase in the number of particles, which reduces the distance between them owing to the limited space available for particle movement.
(ii) Particle fracture: The JH2 damage parameter for the particles is found to be zero for both volume fractions (10% and 20%) in the two types of FGMMCs with the power law and Gaussian distributions.The distribution of the JH2 yield stress www.nature.com/scientificreports/parameter for the particles is shown in Fig. 16, which indicates critically loaded particles that are susceptible to damage.
For samples with a power law distribution (represented in Fig. 16a,b), the range of the yield stress parameter is between 3958 MPa and 4971 MPa , with an average value of 4309 MPa for a VF r of 10% ( n = 2.65 ), as shown in Table 5.For a VF r of 20% ( n = 0.65 ), the range is between 3994 MPa and 5207 MPa , with an average value of 4488 MPa (see Table 5).
For samples with a Gaussian distribution (represented in Fig. 16c,d), the range of the yield stress parameter is between 4104 MPa and 5651 MPa , with an average value of 4530 MPa for a VF r of 10% (see Table 6).For a VF r of 20%, the range is between 4040 MPa and 4993 MPa , with an average value of 4517 MPa .These findings provide  www.nature.com/scientificreports/insight into the distribution of the yield stress parameter, highlighting the potential susceptibility of critically loaded particles to damage in both the power law and Gaussian FGMMCs.
(iii) Interfacial decohesion: For samples with a power law distribution and 10% VF r ( n = 2.65 ), there are three samples ( S 2 , S 5 , and S 6 ) in which complete interfacial decohesion is observed.One of these samples ( S 6 ) is shown in Fig. 17a.Two samples ( S 4 and S 8 ) exhibited interfacial decohesion when the VF r is 20% ( n = 0.65 ), with one of the samples ( S 4 ) shown in Fig. 17b.Regarding the Gaussian distribution, interfacial decohesion is observed in three samples ( S 3 , S 6 , and S 7 ) with a VF r of 10%, and one sample ( S 7 ) with a VF r of 20% .Figure 17c,d shows one sample for each VF r .These findings indicate that interfacial decohesion is a prevalent phenomenon in FGMMCs regardless of the distribu- tion type and VF r .

Comparative analysis of the power law and Gaussian distribution
To compare the two distribution approaches for FGMMCs (the power law and Gaussian distributions), several factors are considered.First, the HRB values are examined, which indicated that the Gaussian distribution approach is the preferred method, particularly for a VF r of 10%.Additionally, to obtain a fair understanding of the preferred distribution approach, the damage incurred by the three different modes the FGMMCs is studied.Because damage does not occur uniformly across all modes, a normalization procedure is conducted for each damage parameter 50 .For each VF r , the damage parameter, or the parameter indicating the occurrence of damage, is normalized based on the results obtained from the samples distributed according to both the power law and Gaussian approaches.The following normalization formula of 50 is used Here, X norm represents the normalized value, X denotes the set of values used for normalization, X min represents the minimum values within X , and X max represents the maximum values within X .By employing this normaliza- tion process, a fair comparison between the damage parameters obtained from the two distribution approaches is achieved, allowing for a complete evaluation of the preferred distribution approach for FGMMCs .Figure 18 presents the normalized values obtained for each distribution approach and VF r .
Based on the analysis of Fig. 18a,c,e for a VF r of 10%, it can be concluded that FGMMCs with a Gaussian distribution are more susceptible to damage, particularly interfacial decohesion (as depicted in Fig. 18e).This is because the higher susceptibility to damage in the Gaussian distribution of samples with a particle VF r of 10% can be attributed to the fact that 60% of these samples exhibited more than 40% damage.In contrast, for samples following a power law distribution, the percentage of samples that experienced damage is lower (approximately www.nature.com/scientificreports/30%).Therefore, for a VF r of 10%, the Gaussian distribution offers higher hardness, but is also more susceptible to damage, although it does not necessarily ensure the occurrence of such damage.
From the examination of Fig. 18b,d,f, it is evident that the Gaussian distribution exhibits lower susceptibility to damage.Hence, for a VF r of 20%, a Gaussian distribution is a favorable choice to achieve both high hard- ness and improved resistance against damage.Considering these results, decisions regarding the distribution approach should consider the specific requirements of the application and evaluate the desired hardness against an acceptable level of potential damage.

Conclusion
A novel finite element model based on a random microstructure is proposed to investigate the deformation and damage behavior of functionally graded metal matrix composites subjected to indentation loading.This study focuses on the influence of the volume fraction and distribution percentage of Al 2 O 3 particles on the thickness of the composite structure.To the best of our knowledge, this study represents the first investigation into the influence of particle volume fraction and distribution on the damage mechanisms of AA6061-T6/Al 2 O 3 func- tionally graded metal matrix composites while considering various damage modes occurring during the indentation process.The model incorporated the Gurson-Tvergaard-Needleman (GTN) damage model to simulate the elastoplastic behavior and damage of the matrix, the JH2 cracking model to represent particle fracture, and the cohesive zone model (CZM) to account for matrix-particle interfacial decohesion.The model was employed in ABAQUS/Explicit software using a Python script.The significant findings of this study are as follows: • The finite element model effectively captures the three primary damage mechanisms in PRMMCs , particularly in FGMMCs , thereby providing valuable insights into optimizing their performance.• The hardness value ( HRB ) increases as the particle volume fraction and size increase.Additionally, the prob- ability of fracture for the Al 2 O 3 particles increases with an increase in the volume fraction.• Al 2 O 3 particles enhance the resistance of FGMMCs to deformation, although their effectiveness is influenced by the particle volume fraction.• PRMMCs can be represented using a variation in the particle concentration throughout the thickness of the structure, resulting in FGMMCs .FGMMCs exhibit a higher hardness value ( HRB ) than traditional Metal Matrix Composites ( MMCs).• This study demonstrates that with a volume fraction of 10% particles and a Gaussian distribution for particle variation through the structure thickness, the HRB is approximately equivalent to the HRB achieved with a volume fraction of 20% particles in MMCs without any variation in particle distribution through the thick- ness.
• FGMMCs with a Gaussian distribution yield a higher HRB than those with a power-law distribution and exhibit closer alignment with the experimental distribution functions.

Figure 1 .
Figure 1.Graded distribution of reinforcing phase in matrix from outer to inner surface based on experimental studies.

Figure 3 .
Figure 3. Flowchart of the present analysis.

Figure 4 .
Figure 4. Finite element mesh for the studied model.

Figure 5 .
Figure 5. Applied force as a smooth step in ABAQUS.

Figure 8 .
Figure 8.The distribution of Al 2 O 3 particles across the FGMMCs structure thickness using power law.

Figure 9 .
Figure 9. Compositional profiles of the FGMMCs using Gaussian distribution.

Figure 10 .
Figure 10.The distribution of Al 2 O 3 particles across the FGMMCs structure thickness using Gaussian distribution.

Figure 11 .
Figure 11.The variations of HRB for different Al 2 O 3 particles VF r and variation through the FGMMCs thickness.

Figure 12 .
Figure 12.Effects of VF r of Al 2 O 3 particles on load-displacement curves and deformed indentation surfaces.

Figure 13 .
Figure 13.Effect of power law distribution on residual von-Mises stress and effective plastic strain during loading.

Figure 14 .
Figure 14.Effect of Gaussian distribution on residual von-Mises stress and effective plastic strain during loading.

Figure 15 .
Figure 15.Void volume fraction in AA6061 − T6 matrix for power law and Gaussian during loading.

Figure 16 .
Figure 16.Yield stress distribution in Al 2 O 3 particles for power law and Gaussian during loading.

Figure 18 .
Figure 18.Comparison between the two approaches of FGMMCs distribution.

Table 1 .
Material constitutive models parameters utilized in this study.

Table 2 .
Parameters employed for CZM of matrix/particle interface.

Table 4 .
HRB results for various particles VF r and variation across the MMCs thickness.

Table 5 .
Output values for FGMMCs using power law distribution.

Table 6 .
Output values for FGMMCs using Gaussian distribution.